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We investigate coupled identical phase oscillators with scale-free distribution of coupling strength. 
It is shown that partially locked states can occur due to the inhomogeneity in coupling and some 
properties of the coupling function. Various quantities of the partially locked states are computed 
through a self-consistency argument and the values show good agreement with simulation results. 
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I. INTRODUCTION 



Synchronization of coupled oscillators is important 
and has been widely studied in varietyof systems from 
physics, chemistry and biology [3, SH, IS [fit- Par- 
tially locked states are states of imperfect synchroniza- 
tion with locked subpopulation and drifting subpopula- 
tion, and have been observed in coupled oscillators with 
distributed intrinsic frequencies [3|, U, l5| . The formation 
of locked subpopulation with the increase of the coupling 
strength is one of the explanations for the transition be- 
tween asynchronous states and synchronous states. In 
these cases, partial locking is due to the inhomogeneity 
in intrinsic frequencies. For some oscillators, the coupling 
strength is strong enough to make the oscillators locked 
to the coupling force by overcoming the desynchroniz- 
ing effect of frequency difference, but for others it is not 
strong enough and they drift. 

However, there can be desynchronizing factors other 
than ones intrinsic to uncoupled oscillators. Recent stud- 
ies of networks of real systems show connection topol- 
ogy and coupling strength can be far from homogeneous 
0,1, BE!- Especially, in many real systems including 
the world-wide web, the Internet, social and biological 
networks, the number of connections per node or the to- 
tal coupling strength per node follows scale-free distri- 
bution (or power-law distribution) P(x) ~ x~ 7 @, [![• 
Theoretical studies have shown that the inhomogeneity 
in coupling can induce asynchronous states or make syn- 
chronization harder to achieve [lj, El, EH, EH, EH, EB| • 

In this paper, we study the dynamics of oscillators 
with coupling strength which has a scale-free distribu- 
tion. It is shown that partially locked states can occur 
due to the cooperation of the inhomogeneity and the cou- 
pling function. In contrast to the previously studied par- 
tially locked states in systems with distributed intrinsic 
frequencies, the partially locked states due to coupling 
inhomogeneity can be bistable with synchronous states. 
Using a self-consistency argument, we analytically obtain 
various values of the partially locked states. 



II. MODEL AND SIMULATIONS 

To focus on the coupling inhomogeneity effect, we con- 
sider systems of coupled identical limit cycle oscillators 
with same type of coupling. In the case of weak coupling, 
the systems can be reduced to the following phase-only 
ones [1,1,1, i]. 



1 N 

J'=l 



0, i = l,2,...,N, (1) 



where 0{(t) is the phase of oscillator i at time t, lo is 
the natural frequency of the oscillators, and N is the 
total number of oscillators. is the coupling strength 
from oscillator j to oscillator i and Kij > 0. H(9) is 
the coupling function obtained by the phase reduction 
method for pair- wise interaction [H, 0, [H, Q . 

Recently, we introduced the following mean-field model 
as an approximation of the model of Eq. dXJ) [l6| . 



K N 

j=i 



i), i = l,2, ...,JV, (2) 
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where Ki (> 0) corresponds to the average coupling 
strength to oscillator i. This model is a simple gener- 
alization of the Kuramoto model [H, 0, @| where Ki is 
the same for all the oscillators, and the coupling inho- 
mogeneity is incorporated in Ki. Due to the mean-field 
coupling, this model, like the Kuramoto model, is easy to 
simulate and analyze. In the following sections, we use 
this model to study the scale-free coupling inhomogene- 
ity and relate the results with those obtained from the 
simulations with scale- free networks @. 

Here, we use H{&) = sin(0 — (3) + Co with Co = sin/3 
and j3 £ K),7r/2) instead of H(9) = sin (2 of the Kuramoto 
model [3, 0, [f|. In most of the previous studies with 
coupled phase oscillators on scale-free networks, models 
with H{6) = sinO have been studied [13, El UM- But 
coupling functions of the form H(6) = sm(6 — [3) + cq 
are more realistic approximations of those obtained from 
coupled limit cycle oscillators [6|, and cq can affect the 
dynamics significantly in the systems with inhomogene- 
ity in the number of inputs or in the coupling strength 
to an oscillators [ll|, QJ] . H(6) = sin(0 + sin j3 is an 
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approximate coupling function for diffusively coupled os- 
cillators such as gap-junction coupled oscillating neurons 
with which the coupling is zero when the two oscillators 
are at the same point in the phase space [20j . 

With the coupling function H(9) = sm(6 — (3) + sin j3, 
j3 G [0, 7r/2), in-phase synchrony is a solution of Eq. ([2|) 
regardless of the coupling strengths, and the state is lo- 
cally stable since H'(0) > (Theorem 3.1 of Ref. (HJ). 

We perform numerical simulations of Eq. ([2]) with 
H(9) = sm(6* — 0) + sin/3 using a fourth order Rungc- 
Kutta method with time step At = 0.01. Unless noted 
otherwise, we use near uniformly incoherent initial con- 
ditions: 8i(0) is chosen randomly from [0,27r) for all i. 
Note that near in-phase synchronous initial conditions 
lead to the locally stable in-phase synchrony. Using the 
rejection method [22j, we randomly select Ki according 
to a truncated scale- free distribution g(K). 



CK-~< for K e [A min , A ma J 
0, otherwise, 



(3) 



where if m i n (max) > and C is a normalization factor that 
satisfies the normalization condition J„°° g(K)dK = 1. 

We obtain C = (7-I) (A mil T 7+1 - A ma3 T 7+1 )- Trunca- 
tion in the distribution is introduced to ensure numerical 
stability but simulations with larger A max show that un- 
truncated scale-free distribution gives only quantitatively 
different results. After the assigning of all the values of 
Ki, we renumber the oscillators according to the ascend- 
ing order of coupling strength to clearly see the dynamics 
dependence on Ki. 

To measure the degree of synchrony, we calculate 



1 N 

N ^ 

3=1 



(4) 



R is the order parameter showing the degree of syn- 
chrony: R is between and 1, with meaning uniform in- 
coherence and 1 in-phase synchrony. In the simulations, 
R approaches a stationary value with small fluctuations 
and the time average of R is calculated. 

Figure [1] shows simulation results of Eq. ||5J). In Figs. 
[IJa) and (b), coupling strength Ki and the distribution of 
Ki are plotted, respectively. K mul — 0.1 and A max = 30 
are used in the simulations. Figures QJc) and (d) are the 
snapshot of the phases (9) and the frequency (9) of oscil- 
lators, respectively. While in-phase synchrony is reached 
from almost all initial conditions for the cases with uni- 
form coupling (Ki = K) [23|, > the system evolves to 
a partially locked state with the coupling inhomogeneity 
for some range of f3 value (Figs. [ljc) and (d)). The os- 
cillators are divided into a phase-locked group giving a 
continuous line in the figures and drifting group giving 
scattered dots. Because the oscillator indices are renum- 
bered according to the order of the coupling strength, the 
locked oscillators and the drifting ones are clearly distin- 
guished. With the coupling strength distribution treated 
in this study, it is observed that oscillators with small Ki 
are locked (Figs. [TJa), (c), and (d)). 
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FIG. 1: Inhomogeneous coupling strength distribution and 
partial locking for the system of Eq. with H(8) = sin(# — 
13) + sin [3. N = 2000 and u> = n. (a) Coupling strength K of 
oscillators. Kmi n = 0.1 and A" max = 30 (Eq. @). The inset 
of (a) shows the entire range of K. (b) Coupling strength 
distribution g(K) for (a). g(K) ~ K' 1 with 7 = 2.0. (c) 
Phase of oscillators at a certain time after the system reaches 
a steady state for the case with j3 = 0.44-7T. In this state, 
oscillators with oscillator index i approximately less than 500 
are locked and others drift. Time averaged order parameter 
R ~ 0.657. (d) Frequency of oscillators for the state of (c). 
Those oscillators which have the same frequency are locked 
ones. The insets of (d) shows the entire range of the frequency. 



This is a new phenomenon. While most of par- 
tially locked states discussed previously were due to 
the inhomogeneity in the natural frequencies of oscil- 
lators 0, 0, @, those of this study are mainly due to 
the coupling inhomogeneity. Without the inhomogene- 
ity, the states cannot exist. The formation of partially 
locked states in this system may be related to the for- 
mation of so-called chimera states in nonlocally cou- 
pled identical oscillators with a similar coupling function 
[25l I2II I27I [28j . In chimera states which can be clas- 
sified as partially locked states, phase-locked oscillators 
spatially coexist with drifting ones. In those systems, 
there is no coupling inhomogeneity and thus cq of the 
coupling function has no significant role. In both of the 
cases, partially locked states occur when j3 is near n/2 
[I]! [2|| [27], HI]. We can see more of the similarity in the 
next section through analysis and simulation. Note also 
that the partially locked states are bistable with in-phase 
synchronous states. The bistability between a chimera 
state and an in-phase synchronous state is a similar one 

[HUE S3, Hi. 

For Co 7^ sin /3, in-phase synchrony is not a solution of 
the system and unless Co is large enough to induce uni- 
form incoherence, we can obtain partially locked states 
with the truncated scale-free distribution or with other 
distributions flql. 
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III. ANALYSIS 

In this section, we analyze the partially locked states 
using a self-consistency argument H, 0, [H . 

Before analyzing partially locked states, we need to 
look into the local stability of uniformly incoherent 
states, because uniformly incoherent states are solutions 
of the model and these states can compete with partially 
locked states. In recent paper [l6j], using population den- 
sity analysis, we showed that for the model of Eq. ([2]) 
with H{9) = co+sin(#— /3), the eigenvalues A = ^i—iv—iuj 
determining the stability of an incoherent state satisfy 
the equations 



2cos/3 
2sin/3 



„ /i 2 + (Kc - v) 2 
00 Kg(K)(Kc - v) 



(i 



/i 2 + (Kc Q - v)- 



dK, 
dK, 



(5) 
(6) 



where g(K ) is the distribution for K. 

For the coupling strength distribution of Eq. these 
equations become 



cos /3 
sin/3 



C 
2 

C 
~2 



K n 



Kr, 



-7+1 



fl 2 + (Kc - v) 2 
K-~* +1 (Kc Q - v) 
H 2 + (Kc - v) 2 



dK, 



dK. 



(7) 

(8) 



In the limit of fi — > + , we can find critical cq above 
which uniformly incoherent state is obtained |16|. If we 
take the conditions co-f^mm — v < and coK max — v > 
without which Eqs. ([7]) and © cannot be satisfied, the 
equations become the following in the limit. 



cos/3 
sin/3 



2^o~ 



C 



-7+1 



lim 

p-»o+ 2 JK 



K-f +1 (Kc - v) 
fj, 2 + (Kc -v) 2 



(9) 



dK. (10) 



Figure [2] shows critical cq (c *) as a function of /3. The 
critical values denoted by solid lines are obtained nu- 
merically from Eqs. © and (fTU]) . When c < c *, in- 
coherent states are unstable. For the case of Eq. @ 
with Co = sin/3 and j(> 1.5) studied in this paper, cq 
is less than cq* as shown in Fig. [2] and thus the inco- 
herent states are unstable. Smaller 7 gives ranges of j3 
values with cq — sin /3 > Co*, but for the simplicity of the 
discussion we restrict the cases to those with 7 > 1.5. 
Since incoherent states are unstable in this situation, we 
get states other than incoherent states even when ini- 
tial conditions are near-uniformly incoherent states. The 
system evolves to a partially locked state or an in-phase 
synchronous state from near-uniformly incoherent states. 

Now, let us analyze the partially locked states. Ku- 
ramoto used self-consistency arguments and derived an 
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FIG. 2: Critical Co (co*) for the stability of an uniformly 
incoherent state as a function of (3 for cases with H(6) = 
co + sin(# — P). Other parameter values are the same as in 
Fig. 1. The curves are obtained numerically from Eqs. Q 
and (|10|) . When c = sin/3, an incoherent state is obtained if 
sin/3 > Co*. For all the values of 7 of this figure, sin/3 < co* 
and thus it implies that incoherent states are unstable when 
co = sin /3. 



equation for the order parameter R to analyze his model 
in which oscillators with distributed frequencies are glob- 
ally coupled The idea is to calculate the order 
parameter by calculating the contribution from locked 
subpopulation and from drifting subpopulation in sta- 
tionary states. Because the calculation of the contribu- 
tion from the subpopulations contains the order param- 
eter, this gives a self-consistent equation for the order 
parameter. A similar argument with space-dependent or- 
der parameter was used by Kuramoto and his colleagues 
to analyze recently observed chimera states in nonlo- 
cally coupled identical oscillators [25|, (26[. In chimera 
states, phase-locked oscillators coexist with drifting ones 
[25l \WL l27l HH . We use the same self-consistency argu- 
ment here to analyze our system. In this case R is the 
same for all the oscillators as in the Kuramoto model, 
but the coupling term is different for each oscillator as 
in chimera states because of the dependency of the term 
on the coupling strength Ki. This inhomogeneity in the 
coupling terms causes the splitting of the population into 
two groups - locked group and drifting group. 

We follow the same argument for the analysis of 
chimera states described in detail in Ref. [28| . The main 
difference is that coupling strength distribution g(K) is 
used instead of coupling kernel and the order parameter 
R is the same for all oscillators in our case. 

We assume the limit of infinitely many oscillators and 
focus only on stationary states. Let denote the fre- 
quency of the population oscillation of Eq. ([3]) after the 
system approaches a stationary state and </> = 9 — fit 
represent the phase of oscillators relative to the average 
oscillation. Then we can rewrite Eq. ([2]) using the order 
parameter Re defined in Eq. ^ as follows. 



= Lu-Q + Ki [i?sin($ 
i=l,2,...,N, 



<f>i- f3) + sin /?]. 



(11) 



where $ = — fit. When the system reaches a stationary 
state, R and <i> do not depend on time. The oscillators 
with KiR > \lu — f2 + Ki sin /3| asymptotically approach a 
stable fixed point <pi* obtained from the following equa- 
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tions. 



Q + Ki sin (3 = KiR sin ( 



* + /?) (12) 



and a stability condition for the fixed point 

cos (cpi* - $ + 0) > 0. (13) 

These oscillators are those which are phase-locked at fre- 
quency fl in the original frame. 

The oscillators with KiR < \u> — ft + Ki sin j3\ drift 
monotonically This subpopulation can be described by 
an invariant probability density p(cj>, K) in the station- 
ary state. The invariant probability density should sat- 
isfy the condition p(<f>,K)v — constant, where v is the 
instantaneous frequency 4>. Therefore, we get 



p(<j>,K) 



(uj ~ n + K sin f3) 2 ~K 2 R 2 
2tt\uj - n + KsinP + KRsin(<P - <\> - (3)\ ' 



(14) 



where the normalization constant is chosen such that 
J^ P (^K)d^=l. 

In the rotating frame of the population oscillation with 
the frequency f2, the order parameter contribution from 
locked subpopulation is calculated as follows. 



dK g(K)e 



i<t>*(K) 



D, 



- lf V* / dK g(K) 
Jd, 

K 2 R 2 - (A + K sin f3) 2 + i(A + K sin (3) 

" ~KR~ ' 

(15) 

where A = lo — fl and T>i is the domain with KR > 
| A + Ksin/3\. We use Eqs. (O and (JTSJ) to calculate 

The order parameter contribution from drifting sub- 
population can be calculated by using population density 
p{t/>,K) of Eq. CEU). 



D d JO 



d(j)dK g{K)p{4>, £T)e # 

g(K) 



dK 
D d KR 



(A + Ksmft) 



sgn (Z{K, A)) y(A + K sin /3) z - K 2 R 2 



(16) 



where ^ d(j> p(4>, K)e 1 ^ is calculated using contour in- 
tegration, sgn(a;) is the sign function and Z(K, A) = 
A + Ksm (3. sgn(Z(K, A)) appears here, because it de- 
termines which pole lies inside the contour. T>d is the 
domain with KR < |A + K sin(3\. 



The order parameter i?e l * in the rotating frame is the 
sum of the contributions from locked subpopulation (Eq. 
(fT5f ) and from drifting subpopulation (Eq. US])). Be- 
cause R and $ are independent of K, we obtain 



R z 



le 



9(K) 
K 



(A + Ksin/3) dK 



9(K) 
K 



K 2 R 2 - (A + K sin (if dK 



9(K)s g n(Z)/ (A + Ks . mp)2 ^ K2R2dK 



(17) 



where T> to t is the total range for K. This gives two inde- 
pendent equations for R and A. 

(R = 1,A = 0) corresponding to the in-phase syn- 
chronous state and (R = 0, A) corresponding to the uni- 
formly incoherent state are solutions of Eq. (fl~7f . But 
note that this fact does not guarantee the stability of 
each state. 

We numerically find the solutions (R < 1,A) which 
correspond to the partially locked states of the system. 
In the cases with coupling strength distribution of Eq. 
©, numerical simulations of the model show that K 
with which oscillators are locked are bounded above, A 
is negative, and (R — sin/3) is negative. Based on the 
simulations and the condition for the locking, we can 
guess that T>i = {K : K min < K < Ki = t^^} 
and T>d = {K : K[ < K < K max }. In addition, 
Z(K,A) = A + Ksinf3 is positive in the domain T>d- 
Using the observation from simulations, we numerically 
obtain R and A from Eq. (JTTJ) . 

We also compute the fraction of drifting oscillators. 



f drift 



N drift 



N 



g(K)dK 



g(K)dK 



K, 

kc 1+1 - K m ^ +1 

K ■ - K ~ 7+: 



(18) 



In Figs. GHa), (b), and (c), we plot the values of i?, 
A, and f drift obtained both from simulations and anal- 
ysis for partially locked states. The analysis (solid lines) 
shows good agreement with simulation results (symbols). 
As mentioned, the in-phase state is stable for the entire 
parameter range of this system (not shown in the fig- 
ures) . We use near-incoherent states as initial conditions 
in simulations and get partially locked states for (3 > (3* 
which depends on 7. For (3 < (3* , the system reaches 
to in-phase synchronous state. Simulations and analysis 
show that (3* becomes larger as 7 increases. This is con- 
sistent with the fact that as 7 increases, the coupling dis- 
tribution becomes that of the case with uniform coupling 
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strength and in-phase synchronous states are asymptoti- 
cally reached from almost all initial conditions in the case 
with uniform coupling strength [23|, [HI]. The partially 
locked states exist for (3 values near tt/2, which is similar 
parameter range for the existence of chimera states in 
the systems with nonlocal coupling [13, • As for the 
chimera states, the transition points to partially locked 
states appear abruptly. This contrasts to the continu- 
ous transition between partially locked states and asyn- 
chronous states in svsteins with distributed frequencies 

sua. 

Note that this system has bistability between an in- 
phase synchronous state and an partially locked state. 

To see how the selection of the state between in-phase 
state and partially locked state depends on the initial 
conditions in the presence of inhomogeneity [29j , we look 
at the asymptotic state as a function of the initial condi- 
tions. The initial values for 6i arc randomly chosen from 
[0, 2nr), where r 6 [0, 1]. The corresponding initial value 



3 d0 



sin(7rr) 



We scan 



of R is given by R = ^ f ' 

the interval [0, 1] for r with step size 0.025 to find r below 
which the system goes to the in-phase synchronous state 
with 0.025 intervals. We denote the r value r*. We aver- 
age the corresponding R values (R* values) over states 
from 10 different configurations. In Fig. [3jd), the filled 
symbols show the R* values. In the simulations, we ob- 
tain partially locked states for R < R* and in-phase syn- 
chronous states for R > R* . No other states are obtained 
in our simulations. Open symbols and solid curves are 
the same ones as in Fig. (Ha). Each dashed curve rep- 
resents the other nontrivial branch of the solution of Eq. 
(JTrJ) . The values of R* are close to but below the dashed 
curves. The states for the dashed curves seem to be un- 
stable and seem to act as basin boundaries for at least 
this type of initial conditions. The fact that the R* val- 
ues are close to 1 shows that the basin of attraction for 
the partially locked states is relatively large compared 
to that of the in-phase synchronous states and thus the 
partially locked states can occur rather naturally in this 
system. Synchronous states are stable even with the in- 
homogeneity but the inhomogeneity makes the size of the 
sync basin smaller [291 ]. 



IV. SIMULATION RESULTS WITH 
TRUNCATED SCALE-FREE NETWORKS 

The results of previous sections are applicable in sys- 
tems of coupled oscillators on networks. 



u> H > An 

t = l,2,...,JV, 



sin 



/3) + sin/3] 



(19) 



where oscillator i is influenced by fc^ neighbors with cou- 
pling strength K according to a coupling topology de- 
scribed by an adjacency matrix A. ki is called the degree 
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FIG. 3: Quantities of partially locked states for different 
values of 7 as a function of (3. (a) Order parameter R. (b) 
A = uj — Q.. (c) Fraction of drifting oscillators. In (a), (b), and 

(c) , the symbols represent the average values of each quan- 
tity which are obtained by first averaging the quantity over 
time and then over simulations with different {Ki} configu- 
rations and initial conditions. The fluctuations in the time 
series of each quantity are small compared to the value of 
each quantity. The curves of (a) and (b) are numerically ob- 
tained from Eq. (|17[) . The curves of (c) are obtained from 
Eq. (|18[) using R and A from Eq. I|17p . The error bars indi- 
cate the standard deviation of the time-averaged quantities. 

(d) Initial conditions that lead to partially locked states. The 
open symbols and the solid curves are for R as in (a). The 
filled symbols denote the value of critical initial order param- 
eter Ro* above which the initial conditions evolve to in-phase 
synchronous states. Each dashed curve represents the other 
nontrivial branch of the solution of Eq. tllT|) for each 7. See 
the text for the initial conditions of (d). Other parameter 
values are the same as in Fig. 1. 



of i. We take the element of adjacency matrix Aij = 1, if 
oscillator j influences oscillator i, and Aij = otherwise. 

We consider the cases of truncated scale-free networks 
Q where the degree distribution of ki, P(k), is given by 



P(k) 



Ck for k £ [fc m i n , fc max ] , 
0, otherwise. 



(20) 



When oscillators are randomly coupled to others and fc m ; n 
is sufficiently large, we can use the following approxima- 
tion for Eq. (fT9| . 



JV 



3=1 



AijHi^Oj — 6 

With this, Eq. (|19[) is approximately equivalent to 

- 0) + sin/3] , 



(21) 



Kb N 
3=1 



(22) 



which is Eq. © with K { being ^ 



We simulate Eq. (fT9")l with networks following a given 
degree distribution P(k) of Eq. (|20p . The networks are 
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FIG. 4: Inhomogeneous degree distribution and partial lock- 
ing for the system of Eq. Jl9j with H(0) = sin(0 - 0) + sin /?. 
io = n and K = 1. (a) Degree distribution P(k) ~ fc~ 7 with 
7 = 2.5. N = 5000, fc min = 200 and fc max = 2000 (Eq. ©). 
(b) Phase of oscillators at a certain time after the system 
reaches a steady state for the case with /3 — 0.46-7T. (c) Fre- 
quency of oscillators for (b). (d) Order parameter R. The 
symbols represent numerically obtained time-averaged order 
parameter for each simulation. The curves denote the values 
obtained from the analysis (Eq. (JTTJl) . 

generated as follows. Using the similar method of Section 
II, we randomly select a positive integer k € [fc m ; n , fcmax] 
and assign it to an oscillator among the TV oscillators as 
the degree of the oscillator. We randomly select k oscilla- 
tors as the neighbors. After the network is generated, the 
oscillators are renumbered according to the ascending or- 
der of the degree. For simplicity of generating networks, 
we use directed networks but bidirectional networks do 
not change the results significantly. 



Figure [4] shows the simulation results of Eq. (| 19|) . Fig- 
ure Ufa) shows the degree distribution of the network. 
Figures E^b) and (c) are the snapshots of the phases and 
the frequency of oscillators, respectively. As in the cases 
with inhomogeneous coupling strength shown previous 
sections, this system also shows similar states. In this 
system, the states have near locked oscillators (Figs. HJb) 
and (c)). We compare the simulation results of this sys- 
tem with those obtained from Eq. (fT7f using the fact Ki 
corresponds to Figure (0])(d) shows good agreement 
between them. 



V. SUMMARY AND CONCLUSIONS 

In summary, we have investigated coupled identical os- 
cillators with scale-free distribution of coupling strength 
and found that partially locked states can occur due to 
the inhomogeneity and the coupling function. Various 
quantities of the partially locked states have been com- 
puted through a self-consistency argument. This study 
contrasts with the previous studies in the fact that the 
coupling inhomogeneity instead of the frequency inhomo- 
geneity is the main cause of partial locking and partially 
locked states can be bistable with synchronous states. 
Our findings may help further understanding of syn- 
chronous behavior on inhomogeneous networks. 
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